### REPLICATION FILE -- TAB1, FIGC1 -- STUDY 2, PERCEPTIONS
### Homola, Rogowski, Sinclair, Torres, Tucker, and Webster
### "Through the Ideology of the Beholder: How Ideology Shapes Perceptions of Partisan Groups"
### Political Science Research and Methods

rm(list=ls())
set.seed(12435)
options(stringsAsFactors=F)
library(foreign); library(plyr); library(RColorBrewer); library(survey)

# log file/sink
sink("HRSTTW_Stereotypes_Tab1FigC1_R.txt")

data<-read.dta("HRSTTW_Stereotypes_Study2.dta")

# Recode answers for analysis
selfatt <- c("ATTS1S27","ATTS2S27","ATTS3S27", "ATTS4S27","ATTS5S27","ATTS6S27", "ATTS7S27",
             "ATTS8S27","ATTS9S27","ATTS10S27")
data[,selfatt[1]] <- ifelse(as.numeric(data[,selfatt[1]])<=1,NA,
                            ifelse(as.numeric(data[,selfatt[1]])<=3,1,
                                   ifelse(as.numeric(data[,selfatt[1]])<=7,0,NA)))
for (i in 2:length(selfatt)){
  data[,selfatt[i]] <- ifelse(as.numeric(data[,selfatt[i]])<=-1,NA,
                              ifelse(as.numeric(data[,selfatt[i]])<=2,1,
                                     ifelse(as.numeric(data[,selfatt[i]])<=6,0,NA)))
}


###
### Table 1: Descriptive Statistics: First-Order Beliefs
###

totdes <- svydesign(~1, weights = ~mar2014wt1, data=data)
propsbypid <- rbind(apply(svytable(~ATTS1S27+PID3_MAXN, design=totdes),2,function(x) round((x/sum(x)*100),1))[2,],
                    apply(svytable(~ATTS2S27+PID3_MAXN, design=totdes),2,function(x) round((x/sum(x)*100),1))[2,],
                    apply(svytable(~ATTS3S27+PID3_MAXN, design=totdes),2,function(x) round((x/sum(x)*100),1))[2,],
                    apply(svytable(~ATTS4S27+PID3_MAXN, design=totdes),2,function(x) round((x/sum(x)*100),1))[2,],
                    apply(svytable(~ATTS5S27+PID3_MAXN, design=totdes),2,function(x) round((x/sum(x)*100),1))[2,],
                    apply(svytable(~ATTS6S27+PID3_MAXN, design=totdes),2,function(x) round((x/sum(x)*100),1))[2,],
                    apply(svytable(~ATTS7S27+PID3_MAXN, design=totdes),2,function(x) round((x/sum(x)*100),1))[2,],
                    apply(svytable(~ATTS8S27+PID3_MAXN, design=totdes),2,function(x) round((x/sum(x)*100),1))[2,],
                    apply(svytable(~ATTS9S27+PID3_MAXN, design=totdes),2,function(x) round((x/sum(x)*100),1))[2,],
                    apply(svytable(~ATTS10S27+PID3_MAXN, design=totdes),2,function(x) round((x/sum(x)*100),1))[2,])

diffsbypid <- abs(propsbypid[,2] - propsbypid[,1]) 
totdes <- svydesign(~1, weights = ~mar2014wt1, data=subset(data,PID3_MAXN!="Independent/Other"))
ttestsls <- rbind(c(svyttest(ATTS1S27~PID3_MAXN, design =totdes)$statistic, svyttest(ATTS1S27~PID3_MAXN, design =totdes)$p.value),
                  c(svyttest(ATTS2S27~PID3_MAXN, design =totdes)$statistic, svyttest(ATTS2S27~PID3_MAXN, design =totdes)$p.value),
                  c(svyttest(ATTS3S27~PID3_MAXN, design =totdes)$statistic, svyttest(ATTS3S27~PID3_MAXN, design =totdes)$p.value),
                  c(svyttest(ATTS4S27~PID3_MAXN, design =totdes)$statistic, svyttest(ATTS4S27~PID3_MAXN, design =totdes)$p.value),
                  c(svyttest(ATTS5S27~PID3_MAXN, design =totdes)$statistic, svyttest(ATTS5S27~PID3_MAXN, design =totdes)$p.value),
                  c(svyttest(ATTS6S27~PID3_MAXN, design =totdes)$statistic, svyttest(ATTS6S27~PID3_MAXN, design =totdes)$p.value),
                  c(svyttest(ATTS7S27~PID3_MAXN, design =totdes)$statistic, svyttest(ATTS7S27~PID3_MAXN, design =totdes)$p.value),
                  c(svyttest(ATTS8S27~PID3_MAXN, design =totdes)$statistic, svyttest(ATTS8S27~PID3_MAXN, design =totdes)$p.value),
                  c(svyttest(ATTS9S27~PID3_MAXN, design =totdes)$statistic, svyttest(ATTS9S27~PID3_MAXN, design =totdes)$p.value),
                  c(svyttest(ATTS10S27~PID3_MAXN, design =totdes)$statistic, svyttest(ATTS10S27~PID3_MAXN, design =totdes)$p.value))
totdes <- svydesign(~1, weights = ~mar2014wt1, data=data)

propsbypid2 <- propsbypid
propsbypid <- propsbypid[,c(1,3,2)]
print("Table 1")
data.frame(cbind(apply(propsbypid,1,function(x) paste(x, collapse=" & ")), rep("&",10),
                 diffsbypid, rep("&",10),
                 round(ttestsls,3)))


###
### Figure C.1: Descriptive Statistics: Second-Order Beliefs
###

propsbypid_2OB <- rbind(t(apply(svytable(~ATTS1S30+PID3_MAXN, design=totdes),2,function(x) round((x/sum(x)*100),2))),
                        t(apply(svytable(~ATTS2S30+PID3_MAXN, design=totdes),2,function(x) round((x/sum(x)*100),2))),
                        t(apply(svytable(~ATTS3S30+PID3_MAXN, design=totdes),2,function(x) round((x/sum(x)*100),2))),
                        t(apply(svytable(~ATTS4S30+PID3_MAXN, design=totdes),2,function(x) round((x/sum(x)*100),2))),
                        t(apply(svytable(~ATTS5S30+PID3_MAXN, design=totdes),2,function(x) round((x/sum(x)*100),2))),
                        t(apply(svytable(~ATTS6S30+PID3_MAXN, design=totdes),2,function(x) round((x/sum(x)*100),2))), # Electric car
                        t(apply(svytable(~ATTS7S30+PID3_MAXN, design=totdes),2,function(x) round((x/sum(x)*100),2))), # Marijuana
                        t(apply(svytable(~ATTS8S30+PID3_MAXN, design=totdes),2,function(x) round((x/sum(x)*100),2))), # Ban soda
                        t(apply(svytable(~ATTS9S30+PID3_MAXN, design=totdes),2,function(x) round((x/sum(x)*100),2))), # More taxes
                        t(apply(svytable(~ATTS10S30+PID3_MAXN, design=totdes),2,function(x) round((x/sum(x)*100),2))), # Healthcare
                        t(apply(svytable(~ATTS11S30+PID3_MAXN, design=totdes),2,function(x) round((x/sum(x)*100),2))), # Firearm
                        t(apply(svytable(~ATTS12S30+PID3_MAXN, design=totdes),2,function(x) round((x/sum(x)*100),2))), # Dinosaurs 
                        t(apply(svytable(~ATTS13S30+PID3_MAXN, design=totdes),2,function(x) round((x/sum(x)*100),2))), # Homosexual
                        t(apply(svytable(~ATTS14S30+PID3_MAXN, design=totdes),2,function(x) round((x/sum(x)*100),2))), # Pledge
                        t(apply(svytable(~ATTS15S30+PID3_MAXN, design=totdes),2,function(x) round((x/sum(x)*100),2))), # Fence
                        t(apply(svytable(~ATTS16S30+PID3_MAXN, design=totdes),2,function(x) round((x/sum(x)*100),2))), 
                        t(apply(svytable(~ATTS17S30+PID3_MAXN, design=totdes),2,function(x) round((x/sum(x)*100),2))),
                        t(apply(svytable(~ATTS18S30+PID3_MAXN, design=totdes),2,function(x) round((x/sum(x)*100),2))),
                        t(apply(svytable(~ATTS19S30+PID3_MAXN, design=totdes),2,function(x) round((x/sum(x)*100),2))),
                        t(apply(svytable(~ATTS20S30+PID3_MAXN, design=totdes),2,function(x) round((x/sum(x)*100),2))))[,-c(1,7)]

RepStype_2OB <- propsbypid_2OB[31:45,]
DemStype_2OB <- propsbypid_2OB[16:30,]

orderRep <- NULL
for (i in seq(3,nrow(DemStype_2OB),3)){
  orderRep <- c(orderRep, i-1, i, i-2)
}

RepStype_2OB <- RepStype_2OB[orderRep,]
DemStype_2OB <- DemStype_2OB[orderRep,]


barPlot <- function(freq, y.val, col="gray", brewer=TRUE, rt=FALSE, custom=FALSE, cscol="gray90",...){
  cum.freq <- 0
  for (i in 1:length(freq)){
    tem.freq <- tail(cum.freq,1) + freq[i]
    cum.freq <- c(cum.freq, tem.freq)
  }
  cum.freq <- as.numeric(cum.freq)
  coords.x <- t(sapply(1:5, function(x) c(cum.freq[x], cum.freq[x+1], cum.freq[x+1], cum.freq[x])))
  botbar <- y.val + 0.25
  topbar <- y.val + 0.75
  coords.y <- matrix(rep(c(botbar, botbar, topbar, topbar),5), ncol=4, byrow = TRUE) 
  if(brewer){cols <- brewer.pal(length(freq), col)}
  else{
    if(!custom){cols <- c(col,paste0(col, 1:4))}
    else{cols <- c(col, cscol, paste0(col,2:4))}
  }
  coord.mat <- cbind(coords.x, coords.y, cols)
  if(rt){return(coord.mat)}
  else{apply(coord.mat,1,function(x) polygon(x[1:4], x[5:8], col = x[9], border=FALSE))}
}

### PLOT FOR DEMOCRATS
mygray <- rgb(col2rgb("gray80")[1], col2rgb("gray80")[2], col2rgb("gray80")[3], alpha = 100, maxColorValue = 255) 
pdf("FigureC1_plots_2ndOB_DescDems.pdf")
layout(matrix(c(1,2), ncol=1, byrow = TRUE), heights = c(0.8,0.2))
par(mar=c(2.5, 2.5, 1.5, 0.5), oma=c(0,1,0,1))
plot(0,0, xlim=c(0,110), ylim=c(1,(nrow(DemStype_2OB)+1)), type="n",
     xaxt="n", yaxt="n", xlab="", ylab="")
cols.vec.dem <- rep(c("BuPu","lightcyan", "Blues"),5)
brewer.vec.dem <- rep(c(TRUE, FALSE, TRUE),5)
custom.vec.dem <- rep(c(FALSE, TRUE, FALSE),5)

for (i in 1:nrow(DemStype_2OB)){
  barPlot(DemStype_2OB[i,], i, cols.vec.dem[i], brewer = brewer.vec.dem[i], custom = custom.vec.dem[i])
}
marks <- c(2.5, 5.5, 8.5, 11.5, 14.5)
rightcat.dem <- rep(c(2,3,1,1,4), each=3)
abline(h=marks+1.5, lty=2)
axis(2, at=marks, labels=c("Electric car", "Marijuana", "Ban soda", "Taxes", "Healthcare"),
     tck=0, cex.axis=0.7, cex.lab=0.7, mgp=c(0.5, 0.3, 0), ylab="Democrat misperceptions")
axis(1, tck = -0.3, lwd = 0, cex.axis = 0.7, cex.lab = 0.9, mgp=c(0.5, 0.3, 0))
title(ylab = "", line = 1.5, cex.lab = 0.8, xlab="% Perceived agreement")

### Highlight the border of the real category bar
coordTest <- NULL
for (i in 1:15){
  temp <-barPlot(DemStype_2OB[i,], i, cols.vec.dem[i], brewer = brewer.vec.dem[i], rt=TRUE)
  temp <- temp[rightcat.dem[i],]
  polygon(as.numeric(temp[1:4]), as.numeric(temp[5:8]), lwd=2)
  coordTest <- rbind(coordTest, temp)
}

names.bar <- rep(c("REP", "IND", "DEM"), 5)
### Shade the most popular region by category and PID group
for (i in 1:15){
  temp <-barPlot(DemStype_2OB[i,], i, cols.vec.dem[i], brewer = brewer.vec.dem[i], rt=TRUE)
  temp <- temp[which(DemStype_2OB[i,]==max(DemStype_2OB[i,])),]
  polygon(as.numeric(temp[1:4]), as.numeric(temp[5:8]), density=15, angle = 45,border = NA, col = "black")
  text(103, as.numeric(temp[6])+0.25, names.bar[i], cex=0.7)
}

par(mar = c(0,0,0,0))
paldems.d <- brewer.pal(5, "Blues")
palreps.d <- brewer.pal(5, "BuPu")
palinds.d <- c("lightcyan", "gray90", paste0("lightcyan", 2:4))
plot(0,0, type="n", axes = FALSE, xlab = "", ylab = "")
legend("center", ncol = 6, 
       legend = c("Democrats:",  "Independents:", "Republicans:", 
                  rep("0-20%",3), rep("21-40%",3), rep("41-60%",3), rep("61-80%",3),rep("81-100%",3)),
       cex=0.7, fill=c("white","white","white", paldems.d[1], palinds.d[1], palreps.d[1], 
                       paldems.d[2], palinds.d[2], palreps.d[2],
                       paldems.d[3], palinds.d[3], palreps.d[3],
                       paldems.d[4], palinds.d[4], palreps.d[4],
                       paldems.d[5], palinds.d[5], palreps.d[5]),
       border=rep("white",18))
dev.off()

### PLOT FOR REPUBLICANS
pdf("FigureC1_plots_2ndOB_DescReps.pdf")
layout(matrix(c(1,2), ncol=1, byrow = TRUE), heights = c(0.8,0.2))
par(mar=c(2.5, 2.5, 1.5, 0.5), oma=c(0,1,0,1))
plot(0,0, xlim=c(0,110), ylim=c(1,(nrow(RepStype_2OB)+1)), type="n",
     xaxt="n", yaxt="n", xlab="", ylab="")
cols.vec.rep <- rep(c("Reds", "snow", "RdPu"),5)
brewer.vec.rep <- rep(c(TRUE, FALSE, TRUE),5)
custom.vec.rep <- rep(c(FALSE, TRUE, FALSE),5)

for (i in 1:nrow(RepStype_2OB)){
  barPlot(RepStype_2OB[i,], i, cols.vec.rep[i], brewer = brewer.vec.rep[i], custom = custom.vec.rep[i])
}
rightcat.rep <- rep(c(3,2,2,4,3), each=3)
abline(h=marks+1.5, lty=2)
axis(2, at=marks, labels=c("Firearm", "Dinosaurs", "Homosexuality", "Pledge", "Fence"),
     tck=0, cex.axis=0.7, cex.lab=0.7, mgp=c(0.5, 0.3, 0), ylab="Republican misperceptions")
axis(1, tck = -0.3, lwd = 0, cex.axis = 0.7, cex.lab = 0.9, mgp=c(0.5, 0.3, 0))
title(ylab = "", line = 1.5, cex.lab = 0.8, xlab="% Perceived agreement")

coordTest <- NULL
for (i in 1:15){
  temp <-barPlot(RepStype_2OB[i,], i, cols.vec.rep[i], brewer = brewer.vec.rep[i], rt=TRUE)
  temp <- temp[rightcat.rep[i],]
  polygon(as.numeric(temp[1:4]), as.numeric(temp[5:8]), lwd=2)
  coordTest <- rbind(coordTest, temp)
}

### Shade the most popular region by category and PID group
for (i in 1:15){
  temp <-barPlot(RepStype_2OB[i,], i, cols.vec.rep[i], brewer = brewer.vec.rep[i], rt=TRUE)
  temp <- temp[which(RepStype_2OB[i,]==max(RepStype_2OB[i,])),]
  polygon(as.numeric(temp[1:4]), as.numeric(temp[5:8]), density=15, angle = 45,border = NA, col = "black")
  text(103, as.numeric(temp[6])+0.25, names.bar[i], cex=0.7)
}

par(mar = c(0,0,0,0))
paldems.r <- brewer.pal(5, "RdPu")
palreps.r <- brewer.pal(5, "Reds")
palinds.r <- c("snow", "gray90", paste0("snow", 2:4))
plot(0,0, type="n", axes = FALSE, xlab = "", ylab = "")
legend("center", ncol = 6, 
       legend = c("Democrats:",  "Independents:", "Republicans:", 
                  rep("0-20%",3), rep("21-40%",3), rep("41-60%",3), rep("61-80%",3),rep("81-100%",3)),
       cex=0.7, fill=c("white","white","white", paldems.r[1], palinds.r[1], palreps.r[1], 
                       paldems.r[2], palinds.r[2], palreps.r[2],
                       paldems.r[3], palinds.r[3], palreps.r[3],
                       paldems.r[4], palinds.r[4], palreps.r[4],
                       paldems.r[5], palinds.r[5], palreps.r[5]),
       border=rep("white", 18))
dev.off()

sink()